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Introduction 



T3. 

This work stems from calculations for Hubbard |l], [| || and Anderson lattice j|, |5| 

models in a self-consistent conserving Green's function scheme []6], 0] known as the the 

fluctuation exchange approximation (FEA) ||. For the 2D Hubbard model, special 

features of band structure, such as Fermi surface nesting M and van Hove singularities 



near the Fermi surface [10], [11]] , lead to anomalous frequency and momentum dependences 
of the self-consistent self-energy [[| |TJ|. At half filling the FEA self-energy develops a 
frequency dependence similar to that proposed for a marginal Fermi liquid ]n| , and the 
spin-fluctuation propagator appears to move exponentially close to an instability with 
decreasing temperature. When the spin-fluctuation propagator is sufficiently close to 
this instability, we have been unable to obtain stable converged solutions. For the half- 
filled 3D Hubbard model, where an antiferromagnetic phase transition is expected at 
finite temperature, we have studied the fully self-consistent spin response to a staggered 
magnetic field [I3|. The results are qualitatively similar to those in 2D jT5[, and show no 



magnetic order for a range of U and T well within the antiferromagnetic phase expected 
from Quantum Monte Carlo simulations [IB|. 



For the Anderson lattice model, we have observed the evolution of a coherent quasi- 
particle state with decreasing temperature along with a substantial enhancement in 
the effective mass H. As in the Hubbard model, for parameters relevant to correlated 
electronic systems, the spin fluctuation propagator is very near an instability at low 



temperatures. We would like to calculate reliably the entropy and specific heat as a 
function of temperature, to elucidate the apparent transformation from a lattice of local 
moments in a sea of 'ordinary' conduction electrons to a band of highly renormalized 
quasiparticles. 

The only numerical approximation (other than finite machine precision) in our pre- 
vious implementations of the FEA has been the treatment of the high-frequency tails of 
Green's functions, self-energies, etc. Although we are primarily interested in understand- 
ing the low-energy excitations, high-energy processes make important contributions to 
effective masses, susceptibilities, total energies, etc., and numerical approximations in 
treating these processes must be controlled, and so far as possible eliminated, to ob- 
tain reliable results for the problems described above. In this paper we describe a new 
approach to handling high-frequency tails. We decompose the single particle Green's 
function into two parts, 

G(k,e n )=g(k,e n )+G(k,e n ), (1) 

where, as we will make more precise below, G contains only 'low' frequency parts of G 
and g contains the leading 'high' frequency parts (through some order). The important 
observation is that relatively little information is contained in the high frequency tails of 
the Green's function. The crucial trick is to find an analytic form for g that describes the 
high frequency behavior of G accurately and leads to tractable analytic expressions for 
the contributions from g to susceptibility bubbles, T-matrices and self-energies. Most of 
the detailed information about correlations resides in G, which is much less sensitive to 
the frequency cutoff than was the original G. In this language, most previous approaches 
to solving the FEA numerically correspond to taking g to be identically zero |L7| . 



Taking advantage of massively parallel computers requires scalable algorithms that 
perform efficiently for a wide range of problem sizes using virtually any number of 
processors. To this end, our algorithm solves the equations of the FEA iteratively, 
making use of discrete Fourier transforms at various stages of the calculation to make 
each step embarrassingly parallel. To motivate our new approach, we first sketch a less 
accurate but more straightforward way to calculate the FEA self-energy 

Standard Implementation of the Fluctuation Exchange Approximation 

Central to propagator renormalized perturbation theory is Dyson's equation, relating 
the renormalized propagator G to the self-energy E, 

G-\k, e n ) = Go x (k, e n ) - S(k, e n ), (2) 

where Go is the Green's function of the non-interacting system. For simplicity we will 
discuss only the simplest paramagnetic Hubbard model, and will include in the self- 
energy only the proper second-order diagram and the contribution from exchanged spin 
fluctuations; including the density and pairing fluctuations is straightforward and in- 
volves no new matters of principle. The fluctuation propagator is constructed from a 
susceptibility bubble which is a convolution product of renormalized propagators, 

T 

X P ft(q, w m ) = -— Y, G ( k + q,e n + uj m )G(k, e n ). (3) 

k.n 



In terms of Xph the fluctuation propagator (T-matrix) is simply 

Ux P h(q,Ur, 



T(q,u m ) = - 



m) 



1 - Ux P h(<i, 



w. 



(4) 



and the self-energy is a convolution of the Green's function with the sum of the T— matrix 
and the susceptibility bubble, 

T 

E(k,e n ) = U 2 — J2 G(q + k,u m + e n )lx P h(<l,u m ) + T(q,Lj m )]. (5) 

q,w m 

To obtain the self-consistent self-energy, one starts with a guess for G (e.g. Go) 
and calculates the self-energy from Eqs.(|3|-[5|). The resulting E is then used in Dyson's 
equation to update the propagator, and this procedure is iterated to some level of ap- 
proximate self-consistency 

The sums in these equations extend over all momenta and all frequencies, but in a 
numerical calculation we can include only a finite number of terms. For the momentum 
sums this is simply equivalent to taking a finite lattice with periodic boundary conditions. 
Truncating the frequency sums admits no simple physical interpretation, however, and 
providing an alternative to truncation is the focus of this paper. From the perspective of 
this section, the most straightforward procedure is to introduce a sharp cutoff, by setting 
the Green's function, self-energy, etc. equal to zero for all \e n \ > e c , which we will call 
the sharp cut-off scheme. This leads to highly artificial behavior of the susceptibility, 
T-matrix, and self-energy with increasing frequency. More important from the point of 
view of instabilities and phase transitions, the self-energies and susceptibilities at low 
frequencies lose high-frequency contributions from Green's functions, susceptibilities, 
and T-matrices. 



Posing the Problem Another Way 

The fluctuation-exchange approximation for the Hubbard model (and related lattice 
models such as the Anderson lattice model) has a special feature that makes it especially 
well-suited for a fine-grained SIMD parallel computer such as the Connection Machine: 
the bare interaction is completely local in space and time, and the approximation does 
not introduce any nonlocal effective interactions. As a result, all equations of the the- 
ory can be solved completely in parallel at each point in either (k, e n ) space (Dyson's 
equation and the T-matrix equations) or (r, r) space (susceptibilities and self-energy), 
without the need to evaluate directly any convolutions. In the (r, r) representation, the 
susceptibility bubble is simply 

X pfc (r ) r) = -G(r,r)G?(-r,-r), (6) 

and the self-energy is 

£(r, r) = U 2 [ Xph (r, r) + T ph (r, r)]G(r, r), (7) 

while the natural representations of the T-matrix and Dyson equation are given by 
Eqs. (0) and (Q) above. From this point of view, a possible approach is to begin from 



Go(e n ) with a sharp high-frequency cutoff as before, and to transform back and forth 
between (k, e n ) and (r, r) using fast Fourier transforms (FFTs); we will call this the En- 
scheme. This yields Green's functions at a discrete set of evenly spaced r-points between 
and /3, but the sharp cutoff causes endpoint ringing near r = and r = (3. Because 
many physical quantities come from the Green's function at precisely these endpoints, 
a better approach is to begin the the calculation from the exact G (t) sampled on a 
uniform mesh of r-points, which, after an FFT, has the effect of introducing a gentle 
high-frequency cutoff in G (e n ). The previous calculations described in the introduction 
use this approach, which we will call the r- scheme 0. 

Contributions from High-Frequency Parts of the Green's Function 

Repeated integration by parts of the Fourier integral for G(r, e n ) shows that the disconti- 
nuities of the Green's function and its derivatives at r = determine the high-frequency 
behavior of the Green's function, 

G(r , e „ ) = =^M + .... + Mp^w + tir; y^q^ii^ (8) 

ie n [ie n y +1 [iEn) p+1 Jo Otp +1 

where 

ac«(,)E^'- T » - yG < r - T > . (9) 

••>-" t=o+ dr p ' " 



drP 



T=0" 



Substituting this expression for the Green's functions and a similar high-frequency ex- 
pansion for £ into Dyson's equation leads to expressions for the discontinuities of the 
renormalized propagator in terms of unrenormalized single-particle energies £k an d dis- 
continuities in S, 

G(k - E " ) ~ ir n + i^r + (Sj5 + ' ' ' (10) 

where S^ is the Hartree-Fock contribution to the self-energy For the Hubbard model 
with nearest-neighbor hopping only, Eq. (0) gives the first two discontinuities as 

AG(r) = -5 r , , AG?(r) = (-ii + Un/2)6 T> o-t6\ r \ tl , (11) 

where \x is the chemical potential, t is the hopping matrix element, n is the (self- 
consistent) density and U is the Hubbard interaction. 

We write G as the sum of an analytic part g(r,e n ) containing the leading high- 
frequency behavior, and a part G(t, e n ) represented numerically up to a maximum fre- 
quency e c , 

G(r,e n ) = G(r,e n ) +g(r,e n ). (12) 

A simple analytic form for g(r,e n ) that includes the discontinuities of Eq. (|Tl| ) is 

g{r, e n ) = -AG(r) Q (e n , x (r)) + AG"(r) Q 1 {e n , an(r)) (13) 

with 

1 1 



Qo(e n , x (r)) = - 



ie n -x (r) ie n + x (r) 
4 



- — for e n -^ oo, (14) 

ZCtx 



Qi(e n , xi(r)) 



2xi (r) 



1 



IE, 



1 



for e r 



oo. 



U£n 



(15) 



a;i(r) ie n -xi(r)). 

The discontinuity and derivative discontinuity of G are included in g independent of the 
choices for x (r) and a?i(r). We choose these parameters by setting G(r, 0) for |r| = 0, 1 
and G'(0, 0) equal to zero; we show below that this choice is optimal when forming the 
second-order self-energy 

In Fig. ([!]) we illustrate this decomposition of the non-interacting Green's function 
in both r and e n space. The solid curve is the full G , the dashed line is G — Q , and and 
the dotted line shows the final numerical part, G = G — Q + fiQi, which is represented 
on a discrete r-mesh and transformed with an FFT. The smoother and smaller G has a 
spectral weight that is effectively confined to low frequencies so that the errors introduced 
by e c are much smaller than those incurred in an FFT of the full Green's function in 
any of the standard cutoff schemes. Analytic terms are Fourier transformed exactly and 
functional forms keep track of contributions to infinite frequency. 
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Figure 1: The non-interacting Green's function Go at r = as a function of 
r (left) and its modulus as a function of e n (right) for a ID Hubbard model 
for T = O.lt and fi = 1.21 The solid curve is Gq(v = 0), the dashed curve is 
G (r = 0) - Q , and the dotted line is G(r = 0) = G (r = 0) - Q + fi Qi- The 
removal of discontinuities at r = in the part of Go represented numerically 
corresponds to removing high-frequency tails in e n -space; Go(e n ) has nearly 
vanished for frequencies greater than a bandwidth W = At. Note that all 
energies are measured in units of the hopping matrix element t. 

As described earlier, the FFT method takes advantage of the relative simplicity of 
expressions like that for the second-order self-energy in (r, r) space, 

S 2 (r, r) = -U 2 G(r, r) G(r, r) G(-r, -r). (16) 

Using our decomposition of G, a part of this expression can be calculated analytically, 

o- 2 (r, r) = -U 2 g(r, r) g(r, r) g(-r, -r) (17) 

since this consists of simple functions of r with analytic Fourier transforms. We choose 
xo(r) and Xi(r) so that <72(r, r) contains the leading discontinuities of S(r, r). This 



choice is optimal in that the remaining numerical piece is continuous to second order at 
t = 0. For example, the leading discontinuity in X 2 (r, r) is given by 



A£ 2 (r) = £ 2 (r,0+)-£ 2 (r,(T) 

= -£/ 2 G(0,(T)G'(0,0 + )AG'(r), 



(18) 



which is finite only for r = 0. Since xq and X\ are chosen such that G(0, 0) vanishes, 
this discontinuity is contained entirely in a 2 (r,r), 



A£ 2 (r) = -U 2 0(0, 0") g(0, 0+) Ag(r) = Aa(r) 



(19) 



In Fig. (H), we show E 2 calculated with G = Gq as a solid curve and the numerical part, 
S 2 — cr 2 , with a dotted curve. It is only this numerical function whose Fourier transform 
is approximated with an FFT. As shown in Fig. (^) this leads to a large reduction in 
error in the self-energy in (k, e n ) space at all frequencies in comparison to the errors 
found with the sharp cut-off, r— , and e n - schemes. 
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Figure 2: Full second-order self-energy (solid) calculated with the bare prop- 
agator Go for a ID Hubbard model with U = t, n — —0.7t and T = 0.04t and 
the the numerical part (dashed) obtained by subtracting an analytic contribu- 
tion evaluated using g. 

In a self-consistent calculation, the optimal parameters are adjusted iteration by 
iteration as correlations change the r — > values of the Green's function and its deriva- 
tives. In Fig. (f|) we show the Eq point of the self-consistent second-order self-energy 
as a function of the number of points kept in the representation for G compared to 
that obtained with the r-scheme. As this figure illustrates, the improvements realized 
in Fig. (Q) are also realized in a self-consistent calculation. In particular, as shown for 
the self-energy at its lowest frequency, Eq = irT, substantially fewer points are required 
for this method to achieve the infinite frequency cut-off limit (obtained with, say, the 
r-scheme) to acceptable accuracy. 

It is possible to achieve results with higher accuracy and more rapid convergence by 
adding more analytic terms in g(r,r) to include higher-order derivative discontinuities 
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Figure 3: The second-order self-energy calculated using Go for the ID Hub- 
bard model with U = t, \x = 1.4, and T = 0.04. The exact result (solid line), 
is compared with this method (o), the sharp cut-off method (□), the r-scheme 
(A), and the e n -scheme (o). For all frequencies this method leads to signif- 
icantly more accurate results than the traditional frequency cut-off schemes. 
Note that every other point is plotted in the left panel for greater clarity. 

at r = 0. The optimal parameter choices for each function are again determined by 
requiring that the leading discontinuities in products such as the second-order self-energy 
are contained entirely in analytic terms. This requirement is satisfied if the x»(r) are 
chosen so that the discontinuities and the values of G(r,r), G'(t,t), etc. at r = are 
contained in analytic terms to the fullest extent possible with the available parameters. 

For the full fluctuation exchange approximation, we introduce analogous asymptotic 
expansions for the particle-hole and particle-particle T-matrices, and again represent 
analytically the leading asymptotic behavior. The discontinuities of the T-matrix are 
determined by the discontinuities of the particle-hole and particle-particle bubbles, which 
are in turn determined by the single-particle Green's functions. The optimal parameters 
for the analytic part of the T-matrix are again chosen so that the leading discontinuities 
in the self-energy terms E(r, r) = T(r, r)G(r, r) are contained in the expressions which 
can be treated analytically, cr(r, r) = i(r, r)g(r, r). This requirement is satisfied if the 
values of T(0, 0), T'(0, 0), etc. are contained entirely in the analytic part of the t-matrix, 
t(r,T). 



Calculation of Thermodynamic Quantities 



Thermodynamic properties obtained from the grand thermodynamic potential calculated 
in a conserving approximation (such as the FEA) are guaranteed to be consistent with 
those obtained from a direct calculation involving the self-consistent Green's function 
and self-energy || . In general, thermodynamic properties and the grand thermodynamic 
potential are particularly sensitive to the high-frequency behavior of the propagator. A 
familiar example of the sensitivity to high-frequency parts of the Green's function is the 
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Figure 4: The self-consistent second-order self-energy at Eq for the ID Hubbard 
model with U — 4t, T — 0.1 and fi = —0.5 as a function of the number of 
frequency points kept in the numerical part of G. The rate of convergence with 
respect to the number of frequency points obtained with this method (filled 
diamonds) is substantially improved with respect to the ordinary r-scheme 
(open circles). The solid line represents the e c — » oo limit of the ordinary 
r-scheme. 

slowly converging frequency sum that results when the density is calculated by tracing 
the Green's function. It is not surprising that thermodynamic properties have proven 
difficult to calculate since changes in temperature produce only small relative changes in 
quantities like the free energy, which may be smaller than a fictitious temperature depen- 
dence introduced by the handling of the high frequency cut-off. In Fig. ([5]) we show that 
the calculation of the entropy can be achieved keeping a modest number of Matsubara 
frequencies when the method presented in this communication is used. Here, the entropy 
S is computed two ways: (1) by numerically evaluating 5* = —dF(T,N)/dT\ N (open 
symbols) where F is the Helmholtz free-energy, and (2) by evaluating S = (E — F)/T 
where E is the total energy (closed symbols). Again, this method (squares) produces 
more accurate results for a given number of frequency points than the r-scheme (as 
shown in the figure, the latter can even yield unphysical values of the entropy when 
a small number of frequency points are kept in frequency sums). Reliable thermody- 
namic calculations based on self-consistent perturbation theories may prove helpful in 
understanding the thermal properties of interacting quantum systems, especially since 
the small system sizes typical of exact methods make accurate calculations of thermo- 
dynamic properties problematic. 
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Figure 5: Self-consistent calculation for the entropy for the 3D Hubbard 
model with 16 3 sites, T = 0.1, n = 0.5, and U = 4 as a function of the number 
of frequency points. The open symbols are obtained using the formula S = 
—dF(T,N)/dT\]y and the closed symbols from the formula S = (E — F)/T. 
The results for this method (squares) converge much more rapidly than those 
from the r-scheme (circles). 
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